city_code<-'bj'
source('R/load_spatial_data.R')
#source('R/predict_rate.R') #optional  slow 
pty<-11    #产品类型
source('R/pp_sale.R')

#城市级别模型建立
cat(city_code,'model building...','\n')
ppi_sale.plot<-ppi_sale%>%
    # mutate(year=floor(ymid/12),
           # month=ifelse(!ymid%%12,12,ymid%%12))%>%
    group_by(ha_code,x,y,year,month)%>%
    summarise(price=mean(saleprice,na.rm=T))
write.csv(file = '/tmp/ppi_sale.csv',ppi_sale.plot)
var.names.except<-c('city_code','ha_code','proptype','year','month',
                    'salecount','name','ha_cl_code',
                    'ha_cl_name','dist_code',
                    'dist_name')
sale<-ppi_sale%>%dplyr::select(-one_of(var.names.except))
sale$yearmonth<-paste(ppi_sale$year,ppi_sale$month,sep = '-')
date<-sort(unique(sale$yearmonth))
len<-length(date)
fit.vec_city<-vector(mode='list',len)
names(fit.vec_city)<-date
formula.lm<-'log(saleprice)~.-yearmonth-ymid-volume_rate-greening_rate'%>%
    as.formula
library(formula.tools)
rvars<-rhs.vars(formula.lm,data=sale)
formula.lm.2<-paste(lhs.vars(formula.lm)%>%attr('term.labels'),
                   paste(rvars,collapse = '+'),sep='~')%>% as.formula
for (i in 1:len){
    sale.new<-sale%>%filter(yearmonth==date[i])
    fit.vec_city[[i]]<-step(lm(formula.lm.2, data = sale.new),trace=F)
}
sapply(fit.vec_city,function(x) summary(x)$adj.r.square)%>%fivenum
# Linear Regression
# Robust Linear Regression
# Support Vector Machine
# Decision Tree
# Random Forests
# Gradient Boosting Machine
# Multivariate Adaptive Regression Splines

data.reg<-sale%>%filter(yearmonth==date[i])
avm.reg.ols<-step(lm(formula.lm.2,data=data.reg),trace=F)
# tune.lm(lm, formula.lm.2, data = data.reg)
summary(avm.reg.ols)
# svm,random forest, regression tree  ,knn , rpart regression
library(e1071) 
avm.reg.svm<-svm(formula.lm.2,data=data.reg)
avm.reg.rpart<-best.rpart(formula.lm.2,data=data.reg,minsplit = c(5,10,15))
library(MASS)
avm.reg.rlm <- rlm(formula.lm.2,data=data.reg)

library(randomForest)
avm.reg.rft <- randomForest(formula.lm.2,data=data.reg, ntree=50)

library(gbm)
avm.reg.gbm <- gbm(formula.lm.2,data=data.reg,n.trees=50,shrinkage=0.005 ,cv.folds=10)
# best.iter <- gbm.perf(avm.reg.gbm,method="cv")
# train_pred <- predict.gbm(GBM_model,trainset,best.iter)
# test_pred <- predict.gbm(GBM_model,testset,best.iter)

library(earth)
avm.reg.ear <- earth(formula.lm.2,data=data.reg)
# 
# fit.vec_city.gwr<-vector(mode='list',len)
# names(fit.vec_city.gwr)<-paste(
#     floor(date/12),
#     sprintf('%02d',ifelse(!date%%12,12,date%%12)),sep='-')
# library(spgwr)
# library(maptools)
# library(formula.tools)
# rvars<-rhs.vars(formula.lm,data=sale.new)
# formula.gwr<-paste(lhs.vars(formula.lm),
#                    paste(rvars,collapse = '+'),sep='~')%>%
#     as.formula
# for (i in 1:len){
#     sale.new<-sale%>%filter(ymid==date[i])%>%na.omit
#     coordinates(sale.new)<-~x+y
#     proj4string(sale.new)<-'+init=epsg:4326'
#     bw.gwr.vl<-gwr.sel(formula =formula.gwr,data=sale.new,longlat = T)
#     fit.vec_city.gwr[[i]]<-gwr(formula.gwr, data = sale.new,
#                                hatmatrix = T,
#                                bandwidth = bw.gwr.vl)
# }
